
################################################################################
###########             Supplementary Analysis (Appendix)            ###########
################################################################################


####################   Influential Cases (Figure A1)   #########################

#Load the complete database (with Honduras)
load("lapop_pela_select_cent.RData"); baseF=base2; rm(base2)

cat("\014") 
library(car)
library(dplyr)
library(lme4)
library(texreg)
library(lattice)
library(influence.ME)
library(kableExtra)
options(scipen=999)

#Main Model (model 3 in Table 2)
cgm0=lmer(conf_parlamnt ~ 
            eval_econpais_Dummy_gm+ izq_gm+ voto_presid_gm+ conf_interpers_gm+
            eduniv_gm+ urbano_gm+ mujer_gm+ edad_gm+ casado_gm + 
            wave+ 
            fampol_gm+ fampol_prom+ reeleccion_gm+ reeleccion_prom+
            e_gdppc_gm + v2x_corr_gm + desempleo_gm+
            (1|pais.ano)+(1|paisF), data=baseF)
screenreg(cgm0)

##Cook's Distance
inf <- influence(cgm0, "paisF")

#cutoff
4/17

#Plot: Cook's Distance calculated based on level-2 coefficients
plot(inf, which="cook", cutoff=0.2352941,ylab="Paises", 
     xlab="Cook's Distance", parameters = c(12,18)) 



#####################  Variance Composition (Table A1) #########################

#Data without Venezuela and Honduras
load("lapop_pela_select_cent_wVH.RData")
baseF=base2

###### Null model, country-year #############
m0a <- lmer(conf_parlamnt~1+(1|pais.ano), data=baseF)
sigma2_ind.a <- sigma(m0a)^2
sigma2_enc.a <- lme4::VarCorr(m0a)$pais.ano[[1]]

#Individiduals variance
var.ind.a<-sigma2_ind.a / (sigma2_ind.a + sigma2_enc.a)# 0.9125809
#Country-year variance
var.enc.a<-sigma2_enc.a / (sigma2_ind.a + sigma2_enc.a) #0.0874191

mod.country_year <-c(fixef(m0a),var.ind.a, var.enc.a, NA,NA, AIC(m0a), 
                     BIC(m0a), logLik(m0a)[1])
comp.var<-as.data.frame(mod.country_year)


######## Null model, country #########
m0b <- lmer(conf_parlamnt~1+(1|paisR), data=baseF)
sigma2_ind.b <- sigma(m0b)^2
sigma2_pais.b <- lme4::VarCorr(m0b)$paisR[[1]]

#Individiduals variance
var.ind.b<-sigma2_ind.b / (sigma2_ind.b + sigma2_pais.b)#0.954205
#Country variance
var.pais.b <- sigma2_pais.b / (sigma2_ind.b + sigma2_pais.b) #0.04579496

comp.var$mod.country <-c(fixef(m0b),var.ind.b, NA, var.pais.b,NA, AIC(m0b), 
                         BIC(m0b), logLik(m0b)[1])


######## Null model, wave (year) #########
m0c <- lmer(conf_parlamnt~1+(1|wave), data=baseF)
sigma2_ind.c <- sigma(m0c)^2
sigma2_ola.c <- lme4::VarCorr(m0c)$wave[[1]]

#Individiduals variance
var.ind.c<-sigma2_ind.c / (sigma2_ind.c + sigma2_ola.c)#0.9567506
#Wave (year) variance
var.ola.c <- sigma2_ola.c / (sigma2_ind.c + sigma2_ola.c)#0.04324942

comp.var$mod.wave <-c(fixef(m0c),var.ind.c,NA,NA,var.ola.c,AIC(m0c), 
                      BIC(m0c), logLik(m0c)[1])


######## Null model, country-year + country #########
m0d <- lmer(conf_parlamnt~1+(1|pais.ano)+(1|paisR), data=baseF)
sigma2_ind.d <- sigma(m0d)^2
sigma2_enc.d <- lme4::VarCorr(m0d)$pais.ano[[1]]#Var country-year
sigma2_pais.d <- lme4::VarCorr(m0d)$paisR[[1]]#Var country

#Individiduals variance
var.ind.d<-sigma2_ind.d / (sigma2_ind.d + sigma2_enc.d + sigma2_pais.d)#0.91263
#Country-year variance
var.enc.d<-sigma2_enc.d / (sigma2_ind.d + sigma2_enc.d + sigma2_pais.d)#0.04967
#Country variance
var.pais.d<-sigma2_pais.d / (sigma2_ind.d + sigma2_enc.d + sigma2_pais.d)#0.0376

comp.var$mod.cy.country <-c(fixef(m0d),var.ind.d, var.enc.d, var.pais.d, NA, 
                            AIC(m0d), BIC(m0d), logLik(m0d)[1])


######## Null model, country-year + wave (year) #########
m0e <- lmer(conf_parlamnt~1+(1|pais.ano)+(1|wave), data=baseF)
sigma2_ind.e <- sigma(m0e)^2
sigma2_enc.e <- lme4::VarCorr(m0e)$pais.ano[[1]]#Var country-year
sigma2_ola.e <- lme4::VarCorr(m0e)$wave[[1]]#Var wave

#Individiduals variance
var.ind.e<-sigma2_ind.e / (sigma2_ind.e + sigma2_enc.e + sigma2_ola.e)#0.9140467
#Country-year variance
var.enc.e<-sigma2_enc.e / (sigma2_ind.e + sigma2_enc.e + sigma2_ola.e)# 0.08418377
#Wave (year) variance
var.ola.e<-sigma2_ola.e / (sigma2_ind.e + sigma2_enc.e + sigma2_ola.e)#0.001769544

comp.var$mod.cy.wave <-c(fixef(m0e),var.ind.e, var.enc.e, NA, var.ola.e, 
                         AIC(m0e), BIC(m0e), logLik(m0e)[1])

comp.var[,c(1:5)]<-round(comp.var[,c(1:5)], 3)

comp.var$labels<-c("Intercept", "Individual", "Country-year", "Country", "Wave", 
                   "AIC", "BIC", "Log Likelihood")

comp.var<-comp.var[,c(6,1:5)]
comp.var <- comp.var %>%
  mutate(across(everything(), ~ replace(.x, is.na(.x), "")))

#Table
kbl(comp.var, booktabs = T, "latex", row.names = FALSE,
    col.names = c("","CY", "C", "Y", "CY + C", "CY + Y"),
    caption = "Composition of variance for multilevel models on 'Political trust' with different random intercepts") %>% 
  pack_rows("Proportion of variance", 2, 5, bold=F, italic = T) %>% 
  kable_styling(full_width = F, latex_options ="HOLD_position", font_size = 10.5) %>% 
  row_spec(5, hline_after=T)%>% 
  column_spec(2:6, width="2cm")%>% column_spec(1, width="3.8cm") %>%
  footnote(general = "CY = Country-Year, C = Country, Y = Year (wave)", 
           footnote_as_chunk = T, title_format = c("italic"))



#########   MLMs for a Two-Item Index of Political Trust (Table A2)   ##########

options(max.print = 99999)
rm(list=ls())
load("lapop_pela_select.RData")
base=lapop_pela_select; rm(lapop_pela_select)

#Remove Venezuela
base2 = filter(base, base$paisR != "Venezuela")

#Select variables
base2=dplyr::select(
  base2, id_unico, wave, paisF, paisR, pais.ano, id_lapop, id_pela, year,
  conf_parlamnt, conf_partidos, confpol_2,,eval_econpais, eval_econpais_Dummy, 
  conf_interpers, voto_presid, eduniv, urbano, mujer, edad, casado, 
  izq, familiar_pol_media_w, recambio_media_w, e_gdppc, 
  v2x_corr, desempleo) %>% na.omit()


### Centering LAPOP variables (individual-level) ###

#Left ideology
base2$izq_gm = base2$izq - mean(base2$izq, na.rm=T)

#Presidential turnout
base2$voto_presid_gm = base2$voto_presid - mean(base2$voto_presid, na.rm=T)

#Social trust
base2$conf_interpers_gm = base2$conf_interpers - mean(base2$conf_interpers, na.rm=T)

#economic assessment of the country (dummy)
base2$eval_econpais_Dummy_gm = base2$eval_econpais_Dummy - mean(base2$eval_econpais_Dummy, na.rm=T)

#Urban
base2$urbano_gm = base2$urbano - mean(base2$urbano, na.rm=T)

#Female
base2$mujer_gm = base2$mujer - mean(base2$mujer, na.rm=T)

#Age
base2$edad_gm = base2$edad - mean(base2$edad, na.rm=T)

#Tertiary education (dummy)
base2$eduniv_gm = base2$eduniv - mean(base2$eduniv, na.rm=T)

#Married
base2$casado_gm = base2$casado - mean(base2$casado, na.rm=T)


### Centering PELA variables  ###

#Legislative Turnover to Reelection
base2$reeleccion <- 1 - base2$recambio_media_w
base2$reeleccion_gm = base2$reeleccion - mean(base2$reeleccion, na.rm=T)
base2$reeleccion_prom<- ave(base2$reeleccion, base2$paisF, FUN = mean)

##Relatives in Politics
base2$fampol_gm = base2$familiar_pol_media_w - mean(base2$familiar_pol_media_w, na.rm=T)
base2$fampol_prom<- ave(base2$familiar_pol_media_w, base2$paisF, FUN = mean)


## Centering V-DEM variables ##

#GDP Percapita      
base2$e_gdppc_gm = base2$e_gdppc - mean(base2$e_gdppc, na.rm=T)

#Political corruption
base2$v2x_corr_gm = base2$v2x_corr - mean(base2$v2x_corr, na.rm=T)

## Centering World Bank variables

base2$desempleo_gm = base2$desempleo - mean(base2$desempleo, na.rm=T)


### Analysis ###

#Remove Honduras
base2 = filter(base2, base2$paisR != "Honduras")

#Index fiability
trust_index <- select(base2, conf_parlamnt, conf_partidos)
alpha(trust_index)

cgm0=lmer(confpol_2 ~ 
            eval_econpais_Dummy_gm+ izq_gm+ voto_presid_gm+ conf_interpers_gm+
            eduniv_gm+ urbano_gm+ mujer_gm+ edad_gm+ casado_gm + 
            wave+ (1|pais.ano)+(1|paisF), data=base2)

cgm1=update(cgm0, . ~ . +   fampol_gm+ fampol_prom)
cgm2=update(cgm0, . ~ . +   reeleccion_gm+ reeleccion_prom)
cgm3=update(cgm0, . ~ . +   fampol_gm+ fampol_prom+ reeleccion_gm+ reeleccion_prom)

#Controlling by contextual variables (V-DEM and WB)  
cgm1_a = update(cgm1, . ~ . + e_gdppc_gm + v2x_corr_gm + desempleo_gm)
cgm2_a = update(cgm2, . ~ . + e_gdppc_gm + v2x_corr_gm + desempleo_gm)
cgm3_a = update(cgm3, . ~ . + e_gdppc_gm + v2x_corr_gm + desempleo_gm)
screenreg(list(cgm1_a,cgm2_a,cgm3_a), stars=c(0.1,0.05,0.01))


